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ABSTRACT 

We present a new analytic approach to the disk-planet interaction that is especially useful for planets 
with eccentricity larger than the disk aspect ratio. We make use of the dynamical friction formula to 
calculate the force exerted on the planet by the disk, and the force is averaged over the period of the 
planet. The resulting migration and eccentricity damping timescale agrees very well with the previous 
works in which the planet eccentricity is moderately larger than the disk aspect ratio. The advantage 
of this approach is that it is possible to apply this formulation to arbitrary large eccentricity. We 
have found that the timescale of the orbital evolution depends largely on the adopted disk model in 
the case of highly eccentric planets. We discuss the possible implication of our results to the theory 
of planet formation. 
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1. INTRODUCTION 

The gravitational interaction between a planet and a protoplanetary disk is one of the main topics of the theory of 
planet formation. A low-mass planet embedded in a protoplanetary disk interacts with the disk by the gravitational 
force, and its orbital elements change as a result of the interaction. The change of the semimajor axis of the planet 
is called (type I) orbital migration. It has rec ently been noted that the direction of migration is sensitive to the disk 
model if the planet is in a circular orbit (e.g.. iPaardekooper and M cllcma 2006; Paa rdekooper et alHteOlOallbl ). 

The observations of extrasolar planets have revealed th at there are a number of planets with high eccentricity, and 
the median eccentricity is ~ 0.3 (jUdrv and Santoi 120071) . One interesti ng question here is whether it is possible to 
have an eccentric planet in a circular disk. Recent numerical simulations (Cress well et aLlfeOOTt iBitsch and Klevll2010( ) 
show that the eccentricity always damps. It is reported that the eccentricity damping and the migration timescales 
do not strongly depend on the physical state of the gas (e.g., radiative or locally isothermal) for a planet with high 
eccentricity, and the timescale becomes longer if the eccentricity becomes larger. 

The linear analysis of the interaction between the disk and a pla.net h as been done by a number of authors (e.g., 
IGoldreich and Tremaindll980t I Artvmowiczl 119931 : iTanaka a nd Ward 2004) for a planet with low eccentricity. In par- 
ticular, IPapaloizo u and Larwood (2000) obtained the eccentricity damping timescale and migration timescale for a 
planet with e p> h, where e is the eccentricity of the planet and h is the disk aspect ratio (h = H/r, where H is 
the disk scale height and r is the disk radius), but their approach is restricted to e < 1. Tana k~a~and Wardl (|2004j ) 
performed three-dimensional modified local linear analyses to calculate the gravitational interaction between a disk 
and a planet in an eccentric orbit. They have calculated the density perturbation at every location of the orbit, and 
therefore, it is possible to calculate the instantaneous force acting on the planet at every position on the orbit. The 
instantaneous force acting on the planet is then averaged over the orbital period to obtain the timescales of the orbital 
evolution. We note that since they use the (modified) local approximations, their results can be applied to the case 
where the eccentricity of the planet is small compared to the disk aspect ratio. 

In this paper, we present an analytical model for the interaction between a low-mass planet in an eccentric orbit and 
the disk. In contrast to the previous approach i n which one uses Fourier decom position and calculates the contributions 
from Lindblad and corotation resonances (e.g.. IPapaloizou and Larw ood 200iJ), we make use of the dynamical friction 
formula to estimate the force acting on the planet at every location on the orbit. The force is then averaged over the 
orbital period to obtain the evolution timescales of the orbital parameters. Using this "real-space" model, it is also 
possible to obtain more intuitive pictures of disk-planet interaction. We note that most of the highly eccentric planets 
detected so far are gas giants, but we consider a low-mass planet in this paper to make the problem simpler. Our 
analytic approach would reveal an underlying physical mechanism of the planet-disk interaction, which, we believe, 
would be useful in the investigation of high-mass planets as well. We also note that in the course of the formation 
processes of the planets, it is important to understand the evolution of the orbital parameters of a low-mass protoplanet 
or the core of the gas giants. In this c ase, our approach based on linear perturbation analyses may be directly applicable. 

Our approach is similar to that of ITanaka an d Ward! (|2004| ) in the sense that we calculate the force acting on the 
planet at every instance of the orbit. The major difference, however, is that we use a very simple model for the 
force acting on the planet. As we shall show later in this paper, our approach is especially useful for the case with 
eccentricity larger than the disk aspect ratio, and we expect that it is possible to use our model for a planet with 
an arb itrary large eccentricity. Therefore, we are in the parameter space that is complementary to ITanaka and Ward! 
(2004). However, it is noted that the use of a simple model also poses the limitation of our approach, and we shall 
discuss the applicability of the model towards the end of the paper. 

The dynamical friction in a gaseous medium itself is also an interesting topic investigated by a number of authors. 
iRephaelia nd Salpcter (1980) investigated the stationary pattern around a gravitating object in a homogeneous medium 
using a linear perturbation analysis, and concluded that the dynamical friction force vanishes if the particle's speed is 
subsonic, while the dynamical friction force varies as v~ 2 in the supersonic case, where v i s the speed of the pa rticle. 
The result in the supersonic case is in agreement with the case of the collisionless system (Chandrasckhar 1943). 

The conclusion of zero dynamical friction force seems rather counterintuitive since the drag force may experience the 
su dden drop when the particle's speed becomes from supersonic to subsonic. This apparent contradiction is resolved 
bv lOstrikerl (|1999h . who performed the time-dependent analysis. She showed that the dynamical friction force depends 
linearly with the speed of the particle when it moves at a subsonic velocity, while the force depends on u -2 for the 
supersonic case. The results of lOstrikeil (1999) o btained by linear perturb ation analyses are in good agreement with 
non-linear numerical simulations (|Sanchez-Salcedo and Branden burg 1991). 

After th e se wo r ks, there are a number of works on this topic with different configurations of the problem (e.g., 
IKim et alJ (120051) : IKim and Kiinl (I2007L I2009I) ). but it seems there has not been a work considering a slab geometry, 
which can be applied to the problem of disk-planet interaction. 

In this paper, we first perform a linear perturbation analyses and derive an analytical formula of the dynamical 
friction force exerted on a particle embedden in a homogeneous slab of gas. We then apply the formula to the problem 
of the disk-planet interaction. We consider the case where the planet is in an eccentric orbit. Although the gas flow 
in a protoplanetary disk is not homogeneous, we show that it is actually possible to obtain reasonable results if the 
planet is in an eccentric orbit. We note that the results of the dynamical friction force may have different astrophysical 
applications from the problem of the disk-planet interaction. 

2 Here, "modified" means that they take the terms up to the second lowest order of H/r. 
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This paper is constructed as follows. In Section [5J we derive the dynamical friction force exerted on a point particle 
embedded in a homogeneous gas slab. In Section [3l we apply the dynamical friction formula obtained in Section [2] to 
the problem of disk-planet interaction. We then discuss some possible applications to the planet formation theory in 
Section HI and Section [5] is for summary. 

2. DYNAMICAL FRICTION IN A GASEOUS SLAB 

In this section, we consider the dynamical interaction between a point particle and a homogeneous, viscous gaseous 
slab in which the particle is embedd e d. W e show that it is possible to obtain a dynamical friction formula whose 
behavior is similar to that of lOstrikerl (|1999f ) in this setup, although there is additional dependence on the viscosity as 
well as the Mach number of the particle. We also show that in the limit of a high Mach number, the resulting formula 
does not depend on the viscosity. 



2.1. Basic Setup 

We consider a homogeneous gaseous slab with surface density So- We take the Cartesian coordinate system (x,y) 
in the plane of the slab, and consider a point particle with mass M traveling in the y-direction at a constant speed Vq. 
The basic equations we consider are the (vertically integrated) hydrodynamic equations with a simplified prescription 
of viscosity, 

dV 

" ■ V • (Ev) = 0, (1) 



dt 

dv _ 1 _ „ 



iA7 2 v - Vip, 



(2) 



where E is the surface density, v is the velocity, P is the (vertically integrated) pressure, v is the viscous coefficient, 
and ip is the gravitational potential of the particle. We use a simple form of viscosity in order to keep the problem 
simple and tractable. We use a simple, isothermal equation of state 



P = c 2 E, 



(3) 



where c is the sound speed. We consider an inertial coordinate system in which the particle is always at the origin so 
the background flow velocity is given by v = VQe y . For the gravitational potential, we make use of the form that is 
commonly incorporated in the investigation of the disk-planet interaction in two-dimensional analyses, 



tl> = - 



GM 



yjx 2 + y 2 + e 2 



(4) 



where e is the softening length, which is a sizable fraction of the thickness of the slab in order to mimic the three- 
dimensional effects. 

This form of the gravitational potential ((4]) is the model for the vertically averaged gravitational potential, and 
therefore, the softening length should be of the order of the vertical scale length of the disk. It is to be noted that, 
with this form of the potential, the gravitational potential close to the point mass (typically, the distance closer than 
e) is underestimated. If the vertical averaging is taken, the potential close to the mass behaves as oc logr, where r is 
the distance from the point mass, in contrast to equation ((3]), which behaves as constant for r <J e. 

However, as we shall show later, the main contribution to the dynamical friction force comes from the region r ~ e, 
and the contribution from the region r < e would affect the results by only some factor. It is also noted that the 
dependence of the force on the physical parameters can be captured with the form of the potential given by equation 
([%]), Therefore, we use this form of the potential as a model of the vertically-averaged potential in this paper. The 
outcome of this approximation will be discussed in detail later in this section. We also note that this form of the 
potential is widely used in the numerical calculations of disk-planet interaction. One benefit of using the simple form 
of the gravitational potential given by equation Q is that it is possible to obtain an analytic expression for the 
dynamical friction. 



2.2. Linear Perturbation Analysis 

We now perform the linear perturbation analysis to calculate the surface density perturbation induced by the 
particle's gravity. We use the subscript zero to indicate the background quantities and we denote all the perturbed 



quantities by 5, e.g., surface density is given by E = 
perturbation. Equations ([1} and ^ now become 



<5E. We retain the terms up to the first order of the 



d_6Y, 
dt E 

d 



d 5E 
dy E 

d 



dt 

d_ 

dt 



vo 



dy 



dy 



9 A 

dx 



Sv x 



9 A 

dy 

, d 8E 
dx Eo 

d <J£ 



0, 

dip 
dx ' 



Sv y = -c — — — 

dy E 



(hp_ 
dy' 



(5) 



(6) 



(7) 
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From these equations, we derive a single second-order differential equation for the surface density perturbation, 

I + 2 « - -V 2 (| + ,o|) a - c 2 V 2 a = V 2 ^ (8) 

where we have defined a = 6E/E0 and V 2 = jd/dx) 2 + (d/dy ) 2 is the Laplacian. 

We now consider a steady state where d/dt = 0. In lOstriken (fi"99 9), it is pointed out that it is necessary to perform a 
timc-dependent analysis in order to correctly obtain the dynamical friction force, especially when the particle's velocity 
is subsonic. This is because the contribution to the force coming from the place very far away from the particle is not 
negligible. However, as we shall show below, the analysis assuming the steady state is adequate in the slab geometry 
we consider in this section since the perturbation induced at a distant place from the particle does not contribute to 
the force. 

In Appendix^ we explicitly show that the force coming from the time-dependent terms falls as t . Here, we briefly 
show this by the order-of-magnitude estimate. If we consider the place far away from the particle, the gravitational 
potential is given by 

, (9) 

where r is the distance from the particle. We expect that this gravitational energy is of the same order of the magnitude 
with the perturbed thermal energy of the gas. In the three-dimensional analysis, we therefore expect that 

do) 

where p is the gas density. In the slab geometry, we expect that 

f~^, an 

Eo c z r 

where E is the surface density. In the three-dimensional analysis, the force SF acting on the particle from the gas shell 
at the distance r with the width Sr is 

rri GM5pr 2 Sr (GMfpSr Sr 

6F ~2 ^ ^ — ■ ( 12 ) 

The force from the shell decays only with r _1 . Therefore, when summed over all the shells, the contribution from the 
distant shell is not negligible. In the case of two-dimensional slab geometry, on the other hand, the force from the 
ring at distance r with width Sr is 

._ GMSErSr (GM) 2 E Sr Sr 

" ~ ^2 ~ ^2 ^ ~~2' \ ^) 

The force from the distant ring decays as r~ 2 and therefore, the contribution from the distant ring does not account 
the total force when summed over all the rings. We note here that the contribution to the dynamical friction force 
far away from the planet decays as r" 1 if we integrate all the rings. This explains why the force decays as t~ x in the 
subsonic case if we consi der the s l ab ge ometry. The rings that contribute to the dynamical friction force reside at 
r ~ ct, as pointed out by lOstriken (|1999f ). 

In this paper, we shall show in detail the derivation of the dynamical friction force exerted on the particle embedded in 
a gaseous slab using the two-dimensional approximation. Before going on to the details, we show that the dependences 
of the force on the physical parameters of the gas can be derived in the two-dimensional approximation. 

In the realistic case of the gaseous slab with the thickness of L z (^ e), the interaction between the gas and the 
particle can be well approximated by the two-dimensional analysis for the scales larger than L z (far region). Inside 
this region (near region), the interaction should be calculated in three-dimension. Integrating equation (|12[) between 
the minimum cut-off scale r m i n and ^L z (7 is a factor of the order of unity) , the contribution to the dynamical friction 
force from the near region, Fn, is given by 

Fn „^U2±). (14) 

The contribution to the force from the far region, Ff, is given by integrating equation (|13j) from ^L z to infinity and, 

( GM) 2 E {GMfp 
2 lL z 7c ; 



Ff ~ ~ (15) 



where we have used E ~ pL z . The integral from the far region does not diverge at r = 00 since there is not enough 
mass to contribute to the force in the slab geometry in the far region, in contrast to the case of the homogeneous 
three-dimensional distribution of the gas. It should be noted that the dependences of the force on physical parameters 



3 This is why Coulomb logarithm is involved in the dynamical friction formula. 
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in Fp and Fn are the same except for the logarithmic factor, which only introduces very weak dependences on L z and 
r m ; n . By choosing the appropriate value of 7, it is possible to have an expression that would reproduce the results 
of three-dimensional calculations from two-dimensional calculations. Therefore, we expect that the two-dimensional 
analyses can be used to understand the fundamental aspects of the dynamical friction force. More detailed discussions 
on the two-dimensional approximation, in relation to the application to the disk-planet interaction, can be found in 
Section EH 

We now derive the steady state solution of equation ©. We make use of the Fourier transform to solve the equations, 
but our goal is to derive the expression of the force exerted on the particle, which is, of course, the quantity in the real 
space. For any perturbed quantity f(x,y), we define the Fourier transform by 

f(x,y) = ^ / dk x dk y f{k x ,k v )e^ x+k «y\ (16) 



and the inverse transform is 



f(k x , k y ) = i- / dxdyf(x, y)e~^+ k yy\ (17) 



The steady state of surface density perturbation can be derived from equation © in the Fourier space as 

k 2 t> 

Q = 

c 2 k 2 — v^k 2 + ivvokyk 2 ' 



(18) 



where k 2 = k 2 + k 2 . 

y 

If we know the profile of the perturbed surface density, the dynamical friction force exerted on the particle is 
calculated by 



F v( x ,y) = J d 2 x8Y 1 (x,y)^-. 



(19) 



We note that the ^-component of the force is zero because of the symmetry. From the solution in the Fourier space, 
the dynamical friction force is given by 

/>oo POO 

F y (x,y) = 2 dk y dk x k y <4>lm(a) , (20) 

JO J — oo 

where the Fourier transform of the gravitational potential if) is given by 

V> = — e ~ ke . (21) 

We calculate the dynamical friction force using equations (jT5J) and ([2"0")) . By changing the integration variables first 
from (k x ,k y ) to (fc, 9) via 

k x = kcas6 (22) 

and 

k v = k sin 6, (23) 
and then from k to u — vv^k/c 2 , we can rewrite the integral as 

F y = ^(GMf T r uIe{M 2 iU 2 )e -ru dUi (24) 

C € ^ JO 

where Mq = vq/c is the Mach number of the background flow and T = 1cx.jM.QV. The function Ig(jp,q) which appears 
in the integral is given by 

P 27T sin 2 # 

Ie(p,q)= . 2ma-u . 2 J 9. (25) 

Jo (1— psm Oy+qsm 

It is actually possible to perform the integration involved in Ig (p, q) analytically. The explicit form of this is found in 
Appendix [SJ 

2.3. Expressions in Some Limits 

It is now possible to perform the integral (|24[) numerically to obtain the dynamical friction force in a gaseous slab 
if we give two dimensionless parameters Mq and T. However, for a moment, we look at some limits to obtain analytic 
expressions in order to investigate how the dynamical friction force depends on these parameters. 



6 Muto, Takeuchi, Ida 

2.3.1. Subsonic Limit 

We first look at the subsonic case where Mq <C 1. In this case, we approximate 1 — Mq sin (9 ~ 1 so the integral (|25p 
becomes 

WW CdB-^-. (26) 



Therefore, in the limit of u — > 0, Ie(M§,u 2 ) ~ 7r while in the limit of it — > oo, Iq(Mq,u 2 ) ~ 2tt/u 2 . Connecting these 
two limits, we approximate the integral by 

<27) 

The dynamical friction force is obtained from equation (I24[) . It is possible to perform the integral analytically and we 
have 

* ^ GM ^ ri-d + ffr)e-^ _ ^ ( ^ 



_-uv , y 

y 2 ec 2 



(28) 



where Ei(x) is the exponential integral defined by 



oo e -t 



Ei(-x) = - / dt . (29) 

J x t 

In the limit of subsonic perturbation, we expect that r>l. We then obtain 

7r S (GM) 2 Mo^ 

F «~4^? cT' (30) 

The dynamical friction force is proportional to the speed of the particle and also the viscosity. In the case of an 
inviscid, steady state with a subsonic particle embedded in the gaseous slab, the dynamical friction force vanishes as 
indicated by the time-dependent analysis shown in Appendix [A] 

2.3.2. Supersonic Limit 

We now consider the supersonic limit where Mq ^> 1. We first consider the behavior of the integral Ig(M 2 ,u 2 ) 
separately in the two limits, u <C Mq and u ^> Mq. 
In the case of u <C M 2 , the dominant contribution to the integral comes from 9 ~ 9q, where 6q is given by 

sin 2 #o = ^. (31) 

Approximating sin# = sin(#o + 89) to the first order of 59, 

sin(6» + 59) - sin 9q + 59 cos 9 (32) 

the integral can be approximated by 

4 r°° 1 

'^•''-Wl/'W^WTW*?' (33) 

where the factor of 4 comes from the fact that there are four positions within < 9 < 2tt where sin 2 9 equals 1/Mq, 
and the range of integration is extended from — oo to oo. The integrand of equation (|33l) is the Lorentzian function 
with width ~ u/Mq. The condition where this width must be very small compared to the unity (or more exactly 27r, 
in order to justify the extension of the range of the integration) leads to the condition u <C Mq. Using the formula 

dx 1 = 7T, (34) 

we obtain 

Ie(M*,u 2 )~ (35) 
uM oy /M$ - 1 

In the case of u > Mq, the width of the Lorentzian is too wide to extend the integration range from (0,2n) to 
(—00,00). In this case, we expect that u 2 sin 2 9 ^> (1 — MQsin 2 #) 2 is satisfied in the most part of 9. We also 
approximate 1 — M 2 sin 2 9 ~ —Mq sin 2 9 so we can rewrite the integral as 

Ie(M 2 ,u 2 ) ~ p , 1 (36) 
Jo Mq sm 9 + u 2 
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Using the formula 



we obtain 



dO- 



1 



2tt 



sin^ 9 



I e {M 2 ,u 2 ) 



Wa 2 + 1 



2tt 



(37) 



(38) 



where we have used of the condition u 3> Mq, In summary, in the supersonic limit, the integral (|25[) can be approxi- 
mated by 



Ie(Mlu 2 ) 



2tt 



uMq^M 2 - 1 
2ir 



(u < MqV-M, 



-= (u> Mq^M 2 -!) 



(39) 



We can now perform the integral (|24j) to give the dynamical friction force. The result is 



ec 2 



1 - exp -TMo^/M 2 - 1 



— Ei — YMq\ Mn — 1 



In the supersonic limit, we expect that TMq ^> 1. In this case the expression simplifies to 



Fi, 



,Sq(GM) s 
ec 2 M 2 



(40) 



(41) 



As expected, we obtain that the force is proportional to Vq 2 . 

We note that the dynamical friction force is independent of the viscosity in the supersonic limit. This indicates that 
the dynamical friction force in the supersonic limit is insensitive to the dissipation process. We conjecture that the 
dynamical friction force is insensitive to the physical state of the slab (whether the slab is isothermal or radiative) in 
the case of the supersonic motion. Similar situation happe ns in the discussion of the Lindblad torque exerted on a 
planet embedded in a disk (M eyer- Vernet and Sic ardv 1983). 



2.4. Dynamical Friction Force 

We now come to the point where we consider all the values of Mq. We integrate equation (|24|) numerically to find 
the dependence of the dynamical friction force on the physical parameters. As we have seen in the above analytic 
discussions, there are two dimensionless parameters, T and Mq, in the problem at hand. In this section, we use the 
"Reynoldes number" , Re, defined by 

Re=- (42) 
v 

instead of T, since the parameter T contains both Mach number and the viscous coefficient. 

Figure [T] shows the dependence of the dynamical friction force normalized by Ep(GM) 2 /c 2 e as a function of Mach 
number and Reynoldes number. In the subsonic regime, we see the expected behavior from equation (|30[) where 
dynamical friction force decreases as we decrease the Mach number and viscosity. For the supersonic case, we also 
see the behavior expected from equation (|4T|) where the force decreases as we increase the velocity but the force only 
weakly depends on the values of viscosity. 

Figure [2] shows the dependence of the dynamical friction force on the Mach number in the case of Re = 100, and 
we also show the formulae in the subsonic and supersonic limit, equations ()30[) and (I41|) . These limiting formulae 
well describe the behaviors of the dynamical friction force. Equation (|41[) is actually a good approximation of the 
dynamical friction force even in the case of Mq ~ 2. 

We briefly comment on the divergence at Mq = 1 . In Figures [1] and [2j we see that the dynamical friction force 
diverges at Mq = 1. This divergence comes from the m atching between the background velocity and the sound speed, 
and was already seen in the previous linear analyses by lOstrikerl ([1999) . It is not possible to avoid this divergence by 
the effect of viscosity. We consider that nonlinear effects are important here. 

2.5. Validity of 2D Approximation 

We have derived the dynamical friction force exerted on a particle embedded in a homogeneous gas slab using the 
two-dimensional approximation. We now discuss the validity of this approximation in the subsonic and supersonic 
cases. 

As we have discussed before, the two-dimensional treatment is based on the averaging of the equations in the vertical 
direction, and the phenomena that occur on the scales larger than the vertical averaging scale can be well approximated 
by the two-dimensional calculations. In our model, the vertical averaging scale is given by e, which is the softening 
scale of the gravitational potential in equation (|4|). 
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In the subsonic case, the dominant contribution to the integral (f!M]) comes from Tu ~ 0(1), which is ek ~ O(l). 
This indicates that in the subsonic case, the perturbation with the scale comparable to the gravitational softening 
length becomes important in determining the dynamical friction force. Therefore, it is indicated that the full three- 
dimensional treatment may be necessary to have a more quantitative results. 

In the case of the supersonic motion, let us consider the case where TMq = 2v e/v 1 for simplicity. In this case, 
the integral (|24|) for the dynamical friction force is approximated by 

roo — Tu roo 

F y (x du—j- oc / dke~ ek , (43) 
Jo M o Jo 

and therefore, all the scales satisfying Tu ~ ek < 1 contribute to the force. In other words, all the scales from the large 
scale (small k), where we expect that the two-dimensional approximation is valid, down to the cutoff scale contribute 
to the force. Since the dependence of the cut-off scale e is also present in the case of the supersonic motion, rigorous 
three-dimensional treatment is necessary to give more quantitative expressions of the dynamical friction force. 

However, we expect that it is possible to capture the dependences on physical parameters even in the two-dimensional 
approximation. In the beginning of Section 12721 we have estimated the order of the magnitude of the dynamical friction 
force in equations (|12p and (flU)) . and the contribution from the near region (fT4")) and that from the far region (|15[) 
differ by only a logarithmic factor log(L z /r m j n ). The dynamical friction force derived in this paper corresponds to the 
contribution from the far region. The lack of the logarithmic factor is due to the fact that the potential is smoothed 
in the region r <; e (see also the discussion after equation (j4j ) . We note the similarity of the expressions derived by 
the simple order-of-magnitude estimate (|15|) and more rigorous calculations (jUJ- 

In order to estimate the contribution from the near region, it is necessary to determine the appropriate value of 
r m j n . However, it is not straightforward and t here are several pub lications on this using the non-linear calculations 
(e.g.. IsTanchez-Salcedo and~F jrandcnburg 1999| iKim an d Kim 2009). Let us consider, for a moment, the case of the 
planet embedded in a protoplanetary disk. Naive estimate of r m - m is the radius of the planet, and in the case of an 
Earth-mass planet, it is of the order of 10 4 km. If we use th e typical disk scale h eight H ~ 0.05AU, the value of 
\og{H l VmhO amounts to ~ 7. Recent non-linear calculations (Kim and Kim 2009) suggest that the deviation from 
the linear value may be smaller due to the shock formation. If we use r m , n ~ GM/v 2 , which is the typical scale of 
the distance between the particle and the shock in supersonic motion (e.g.. IKim" a nd Kim 2009 1 ), instead of the planet 



radius for r m ; n , the logarithmic factor is ~ 3 — 6 for an Earth-mass body. We shall give further discussions on the 
non- linear effects on dynamical friction at the end of Section G3 Q In any case, the force derived in this section may be 
different by some factor due to the two-dimensional approximations, but such logarithmic factor can be approximately 
taken into account by choosing the appropriate value for e (or 7 in equation (|15t ). 

Despite the uncertainty of the two-dimensional approximation, especially for the value of e that should be taken, we 
still apply this results to the problem of the disk-planet interaction. It is because many of the calculations to date are 
done in two-dimensional approximations and they involve the potential of the form given in equation (j4]) . One goal of 
this paper is to investigate how well the simple model using the dynamical friction force may be able to describe the 
numerical results of the disk-planet interaction. More rigorous treatment of the dynamical friction force including the 
three-dimensional effects and the vertical stratification will be discussed in future publications. It is also necessary to 
have three-dimensional non-linear calculations of the disk-planet interaction to compare the model. We note that the 
discussion on the form of the potential applies not only to the linear perturbation analyses, but also to the non- linear 
calculations. 

3. GRAVITATIONAL INTERACTION BETWEEN A DISK AND AN ECCENTRIC PLANET 

In this section, we apply the results of the dynamical friction obtained in the previous section to the problem of 
disk-planet interaction. In this paper, we particularly focus on the planet in a highly eccentric orbit. 

Before going on to the main topic, we briefly summarize the notation and the terminology. In this paper, we focus 
on the evolution of the semimajor axis a and the eccentricity e of the planet. We denote the timescale of the evolution 
of the semimajor axis by t a , which is defined by 

(44) 



da I dt 



where da/dt is the time derivative of the osculating element averaged over one orbital period. Similarly, we denote the 
evolution timescale of the eccentricity by t e , which is defined by 

t e = 6 ■ (45) 

\de/dt\ 

In iPapaloizou and Larwoodl ()2000[ ). they use "migration timescale" t m , which is defined by 

*— T77S' (46) 

4 The use of this value for r m ; n instead of the radius of the body is also motivated by the results of the collisionless system, where 90° 
deflection angle appears as r m ; n . 
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where J is the specific angular momentum of the planet. It is noted that t a and t m are not the same, since the 
angular momentum J is given by 

J= v/a(l-e 2 ). (47) 

However, it is possible to express t m in terms of t a and t e . In this paper, we mainly use t a and t e , which we refer to 
as "semimajor axis evolution timescale" and "eccentricity damping timescale" , respectively. We use t m from time to 
time when necessary, and referred to it as "migration timescale" , but it should not be confused with the semimajor 
axis evolution timescale. 

The gravitational interaction between a planet and a protoplanetary disk has been investigated by many au- 
thors. The standard formulation for the linear perturbation analysis of the disk-plane t interaction is done as follows 
(jGoldreich and Tre mainc 1979, 1980; Artymowicz 1993: IPapaloizou and L arwood 2000). The planetary orbit is decom- 
posed into the power series of eccentricity, and then the planetary potential is decomposed into the Fourier series in the 
azimuthal direction. For each Fourier component, the perturbation excited by the planetary potential at resonances 
(Lindblad resonances and corotation resonances) are calculated. Then, the contribution from all the resonances are 
summed up to obtain the force exerted on the planet, which is readily applied to obtain the orbital evolution of the 
p lanet. 

IPapaloizou and Larwoodl (|2000D calculated the torque and energy exchange between a planet and a disk using this 
formalism and obtained the migration timescale and the eccentricity damping timescale. They have found that the 
eccentricity always damps. They have obtained the migration rate and the eccentricity damping timescale as 



t rn = 3.5 x 10 5 /, 1 - 75 

and 

t e = 2.5 x 10 3 / 2 - 



1 + (er /l.SHf 



1 - (er /l.m)4 



1 

1 + T 



ff/roV /2Mj\ /M s \ ( r 



4 \H/r 



0.07 J \ M GD J \M p J \ 1AU 



fH/r \ 4 (2M 3 \ (Mq\( r 



\ 0.07/ \M GD ) \M P J VlAU 



yr (48) 



yr. (49) 



Here, H is the scale height of the disk, ro is the semimajor axis of the planet, Mqb is the disk mass contained in 5AU, 
M p is the planet mass. The parameter f s is related to the softening length e of the planet's gravitational potential by 
f s = (2.5e/H). They performed the analysis for the disk with the constant aspect ratio, H/r, and the surface density 
variation with r~ 3 / 2 . From these equations, the eccentricity damps exponentially when e -C 1, while the damping 
timescale is proportional to e 3 for e > H/r 

The three-dimensional linear perturbation analysis for a planet in an eccentric orbit is done by iTanaka and Wardl 
(|2004f) . They have found, for a small eccentricity e <C 1, the eccentricity damps exponentially as 



te = 1.282-/— 1 - fi-\ (50) 
Af p Erg \r J p 

where M* is the mass of the central star and f2 p is the angular frequency of the planet. If we adopt S = 6 x 
10 2 (r/lAU) -3 / 2 g/cm 2 , which corresponds to the model where the mass contained within 5AU is equal to 2Mj, 
H/ro = 0.07, and M p = M®, equation (JSOJ) gives t e ~ 2.4 x 10 4 yr at tq = 1AU. This value seems one order of 
magnitude larger than equation (|49p with f s = 1. However, given the uncertainty of the softening length, it may be 
possible to obtain the consistent results if one adopts e ~ H . It is necessary to do a detailed comparison between 2D 
and 3D calculations to derive the reasonable values of the softening length. 

The disk-pla net interaction when the pla n et was in an eccentric orbit was investigated by using non-linear num erical 
simulations by Crcsswcll an d Nelsonl (2006); Crcss well et al.1 (|2007| ) and more recently bv lBitsch and Klevl (j2010f ) using 
a fully radiative code. They obser ved that when the eccent ricity was smaller than 0.1, the eccentricity would damp 
exponentially, and the formula by ITanaka and Wardl (|2004D was in good agreement with the numerical results. For 
higher eccentricity such as e ~ 0.3 — 0.4, they observed that the eccentricity damping would slow down, and the 
timesc ale was in excellent agreement with the formula presented bv IPapaloizou and Larwo od (|2000h . iBitsch and Kleyl 
(|2010f) found that for a planet with high eccentricity (e ~ 0.4), full y radiative calculations and locally isothermal 
calculations woul d give the same results. For the migration timescale, iCresswell and Nelson (2006) reported that the 
formula given by IPapaloizou and Larwoodl ([2000) consistently predicted the migration timescale three times shorter 
for the planets with small eccentricity, while for those with large eccentricity, the formula consistently predicted the 
timescale 1.5 times faster. 

So far, the agreement between the numerical calculations and linear analysis is good. However, since this formulation 
involves the expansion of the planetary orbit in a p ower series of eccentricity e, thes e formulae are applicable to the 
case where e <C 1. It is noted that the formula bv IPapaloizou and Larwoodl ((2000) is applicable to the case where 
e ;> H/ro but e <C 1 is still necessary. 

In this paper, we present an alternative model for the disk-planet interaction with a planet in an eccentric orbit, 
which is especially useful for a planet with high eccentricity. We model the interaction using the dynamical friction 
formula and our calculation proceeds as follows. We first calculate the relative velocity between the gas and the planet. 

5 It is noted that t a and t m are different by a factor of two even in the case of the circular orbit. 
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Then, we make use of the dynamical friction formula to obtain the force exerted on the planet at each location of 
the orbit. From these forces, we obtain the evolution of the orbital semimajor axis and eccentricity using Gauss's 
equations, and when averaged over the orbital period of the planet, we finally obtain the timescale for the evolution 
of the orbital parameters. In the following, we describe each step one by one. We shall show that the timescales of 
the evolution of the orbital parameters obtained in this way are in good agreement with the previous work, although 
the instantaneous force exerted on the planet on each location of the orbit may be rather over-simplified. We show 
how such timescales vary as we consider various disk models. We note tha t the u sefulness of the dynamical friction 
formula in the problem of disk-planet interaction was hinted in iPapaloizoul (|2002f l . However, the results were shown 
for some limited number of disk models. 

The big assumption in this model is that we neglect the effects of local shear at the location of the planet. Later in 
this section, we discuss the applicability of this formulation in conjunction with this assumption. 

3.1. Setup of the Problem 

We consider a protoplanetary disk with a central star with mass M* and surface density profile with r -p , 

E = £oW \ (51) 

where So is the surface density at r — r$. We assume for simplicity that the disk is locally isothermal with the 
temperature profile T oc r~ q . Then, the sound speed c of the disk gas varies as c oc r~ q l 2 . We assume that the disk is 
rotating at the Kepler angular frequency and neglect the small difference arising from the pressure gradient. We define 
the scale height of the disk H by H = c/Qk, where c is the sound speed and fix is the Kepler angular frequency. The 
scale height H is therefore varies as if oc r^ 3 " q ^ 2 . The disk aspect ratio, H/r, is given by 

H / r \ (l-9)/2 



= ho - , (52) 

r \r J 

where h$ is the disk aspect ratio at r — r . We consider a planet with mass M p with semimajor axis a and eccentricity 
e. 

In the later sections, we shall often refer to the "fiducial model", which we define as p — 3/2, q = 1, A/* = M Q , 
M p = M®, ho = 0.05, and £o is chosen in such a way that the mass contained within 5AU is 2Mj. In this model, 
the disk aspect ratio is consta nt throughout the disk . We call this model "fiducial" because it is the model used 
in the numerical simulations by ICresswell and Nelson! (|2006f) , which we shall compare our analytic results with their 
numerical calculations. Note that this is not exactly the same as Minimum Mass Solar Nebula. 

The gravitational potential of the planet is given by equation (j4|), where [x, y) is now the coor dinate centered on the 
planet . The softening length e is given by e = eo-ff(r) where eo is the dimensionless parameter. ICresswell and Nelsonl 
(2006) used e = 0.5. We note that in our model, the softening parameter varies as the location of the disk if H varies 
as r. 

3.2. Relative Velocity between the Gas and the Planet 

We first calculate the relative velocity between the planet and the disk. It is readily calculated if we give the planet's 
semimajor axis and eccentricity, since we assume the disk is in Keplerian rotation. 

We assume that the planet mass M p is negligible compared to the mass of the central star M* and use the cylindrical 
coordinate system (r, <fi, z) with the origin at the central star. We denote the mean motion of the planet by n — 
(GM./a 3 ) 1 / 2 , and the true anomaly by /. The velocity of the planet is then 

aen an 

v p = sin/e r H (e cos / + 1) e,*, (53) 

V V 

where e r and are the unit vectors in r- and (^-directions, respectively. The velocity of the gas is 



GM* 

~r~ e ^V^~ e0 - (54) 

The relative velocity between the gas and the planet is calculated by Av = v ff — v p . Now, we define the unit vectors 

(e-r,ejv), where is directed towards the velocity of the planet and is perpendicular to and in the orbital 
plane (see Figure[3|). Unit vectors (e r ,e^) and (e^e^y) are related by 

e r = ex cos /3 — ejv sin /3 (55) 

= ex sin/3 + e^r cos/3, (56) 

where 

cos/3 = - esin / _ (57 ) 
y/1 + 2ecos/ + e 2 
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and 

1 + ecos/ 

sm0 = —== (58) 
sj\ + 2ecos/ + e 2 

In Figure HI we show the evolution of the relative velocity vector Av over one orbital period with e = 0.5. The 
planet feels headwind at the perihelion, and tailwind at the aphelion. Figure [5] shows the amplitude of the relative 
velocity over one orbital period. It is noted that if the planet's eccentricity exceeds the disk aspect ratio, the relative 
velocity is always supersonic, regardless of the position in the orbit. 

3.3. Orbital Evolution of the Planet: Calculation Methods 

We make use of the Gauss's equations to calculate the orbital evolution of the planet. For the evolution of semimajor 
axis and eccentricity, Gauss's equations read 

~dt = an 2 ^ 

and 



de_l 
dt v 



2(e + cos/)T- -sin/iV , (60) 



where v is the speed of the planet and r is the distance between the planet and the central star. We denote the 
perturbing force per unit mass by (T, TV), where T is the component of the force (per unit mass) in the tangential 
direction of the orbit (positive direction is the direction of the velocity), and N is the component in the plane of the 
orbit and perpendicular to T (positive direction is directed inside the orbital ellipse). The evolution equation for the 
true anomaly is given by 



df a 2 nn n \ „ ( ail — e 2 ) „\ , , 

4t = ^T + — Rcosf-S[l + ^ >-smf) , 61 

dt nae |_ V r / 

where rj = (1 — e 2 ) 1 / 2 , R is the component of the force in the radial direction, and S is the component perpendicular 
to R. Components (T, W) and (R, S) are related by 

R = Tcos/3- iV sin/3 (62) 

and 

S = Tsinf3 + Ncosf3. (63) 

For the expressions of the perturbing force (T, W), we use the dynamical friction formula. In this paper, we especially 
focus on the planet with high eccentricity, and therefore assume that the relative speed between the gas and the planet 
is always supersonic. Therefore, we use the limiting form of highly supersonic case, equation (|41[) . The amplitude of 
the force is 

F-&&±, (64) 

where Av is the relative speed between the gas and the planet. The direction of the force is opposite to the relative 
velocity vector. 

In order to find the long-term evolution of orbital parameters, we average equations (|59p . (|60|) . and (|6ip over one 
orbital period assuming the planet is in a fixed orbit. Since the variation of a and e is very small over one orbital 
period, the assumption of a fixed orbit is justified. The averaging is done in such a way that 

~da~ 1 f T ° rb , da 1 f 2rr ,da/dt 



dt T or b Jq dt T or b Jo df /dt ' 

where T or b is the time taken over one period. The averaging is done in the same way for eccentricity. 

3.4. Orbital Evolution of the Planet: Fiducial Model 

In this section, we compare our results with previous numerical calculations for the fiducial model. 
We first show the torque and power exerted on the planet by the disk over one orbit. The power V exerted on the 
planet is given by 

V = v p • F disk = vT (66) 

and the torque T is 

T = e z ■ (r p x F disk ) = rS, (67) 

where r p , v p , and Fdisk are the position vector of the planet from the disk, the velocity of the planet, and the force 
exerted on the planet, respectively. In Figure[6j we show V and T normalized by the planet's energy (the sign inverted) 

E P = (68) 
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and the angular momentum 

L v = VGAf*a(l~e 2 ), (69) 

respectively. 

In the vicinity of the perihelion, the torque and the power are negative because of the headwind towards the planet, 
and they are p ositive in the vic i nity o f the aphelion owing to the tailwind. Similar behavior is observed in the numerical 
simulations of Crcsswcll ct al. (2007). Especially, the behavior of the torque exerted on the planet looks very similar 
except for the small difference between the location of the per ihelion or aphelion and the location of the maximum or 
minimum of the torque (see Figure 8 of ICresswell et all ((20 07)). In our model, this small difference in phase does not 
appear since we assume that the force is proportional to the relative velocity at the location of the planet. 

The comparison of our results with the Figure 10 of ICresswell et all (|2007T ) shows that the behavior of the power is 
very different. This difference can not be attributed to the different sets of parameters they have used from our fiducial 
model. This may be because our model for the force exerted on the planet is very simplified. This difference in power 
leads to the different behavior in da/dt and de/dt over one orbital period, which are shown in Fi gure [71 We note here 
that the evolution of eccentricity is qualitatively different from the Figure 9 of Crcsswcll ct al. (2007). Fortunately, 
however, it is possible to obtain quantitatively the same results with previous numerical calculations when we take the 
average over one orbital period, as shown below. 

Figure [5] shows the results of migration rate t m and the eccentricity damping rate t e . which are obtained by the aver 



aging over one orbital period. Also plotted in this figure are the formula obtained by Papaloizou and Larwooc 



2000) 



given in equations (|48]) and (j49j). For the migration rate, we also show t m given by IPapaloizou and Larwoodl (120000 
times 3/2 for e < 0.5, which describes the results of nu merical simulations better as~no ted by ICresswell and Nelson! 
(2006) . For the damping of eccentricity, the formula by IPapaloizou and Larwoodl (120001) fits the results of numerica l 
simulations well. We see that our results on t e agree very well with the formula by Papaloi zou and Larwoodl ()2000f ) 
for e <J 0.7. For the migration rate, our formulation consistently results in a factor of two slower timescale compared 
to the numerical simulations. Our results deviate from those of IPapaloizou and Larwoodl ([2000) for large values of 
eccentricity. We expect that our treatment is actually better for highly eccentric c ases (discussed in later sections) , 
and than this difference is attributed to the breakdown of the expansion of e used by IPapaloizou and Larwoodl (2000). 

3.5. Orbital Evolution of the Planet: Varying Disk Model 

We have seen that our formulation results in the consistent timescale of the evolution of orbital parameters for the 
fiducial model. We now look at how the timescales behave as we vary the disk model. 
In our framework, the force exerted on the planet can be written in the form, 

Air 

and 

N = Kr-^, (71) 



where K and a are constants. For the model of the force we use in this paper, which is equation (1641) . if we use the 
softening parameter proportional to the disk scale height, e = eoH(r), the constant K is given by 

K = ^M^ > (72) 

and a is given by 

a=p+^-, (73) 

where Ho is the disk scale height at r = rQ. Therefore, models with the same K and a yield the same results. In 
this section, we show the results on how t a and t e depend on the disk profile p, or in more general terms, a. We fix 
the value of surface density to be Eo/(M*/7q) = 10~ 4 , planet mass to be M p /M* — 10~ 6 , disk aspect ratio to be 
Ho/ r o = 0.05, and the softening parameter to be eo = 0.5. The dependence of these parameters on the timescales t a 
and t e is t cx M~ 1 Eq 1 eoH as expected from the form of K, as long as the force exerted by the disk is sufficiently 
small compared to that from the central star. We use the disk model with q = 1, which gives a = p + 1, but we note 
again that the value of a is the key parameter which determines the values of t a and t e . 

Figure [9] shows t a and t e for a = 1, 2, 3,4 (p = 0, 1, 2, 3, respectively for q = 1). In Figure[9l we fix a/r — 1 and see 
how the timescale behaves as we vary the eccentricity. We see that the behaviors of t a and t e strongly depend on the 
imposed disk model. The timescale of the evolution of semimajor axis increases steadily as we increase eccentricity 
if a < 2. If, on the other hand, a > 2, the values of t a first increase as we increase e, but decrease at large e. The 
timescale of the evolution of eccentricity always increases as we increase eccentricity if a < 3, but there is a maximum 
of t e at a certain value of e if a > 3. 

As we decrease the values of a, another interesting behavior occurs for t a . In Figure 1101 we show the behavior of 
t a and t e for a = — 1, 0, 1 (p = —2, —1, 0). We see that the timescale of the eccentricity evolution does not vary much 



Disk-Planet Interaction and Dynamical Friction 



13 



within this parameter range, while t a is negative for almost all the values of e for a = and — 1. The negative sign of 
t a indicates that the direction of the semimajor axis evolution is outward. 

The reason why we obtain the outward semimajor axis evolution can be qualitatively explained as follows. If the 
index p is negative, the surface density increases as a function of radius. As we have seen before, the planet feels the 
tailwind at the aphelion and therefore is exerted positive torque by the gas, while it feels headwind at the perihelion 
and is exerted negative torque. Since the surface density is larger at the aphelion than at the perihelion when p is 
negative, the planet feels overall positive torque, and therefore the semimajor axis increases. 

The disk with negative p seems rather unrealistic. However, the surface density may increase as a function of radius 
locally. The inner edge of the disk is one example where such local increase of the surface density is expected. If a 
planet with a finite eccentricity is placed in such a place , we e xpect that the semimajor axis can increase as a result 
of the disk-planet interaction. Recently, lOgihara et al.l (|2010( 1 suggested that the planets trapped in mean motion 
resonances can be stopped at the disk inner edge. This is because at the gap edge, the planet feels the disk only at the 
place close to the aphelion and therefore it is exerted the positive torque by the disk. This positive torque is balanced 
by the negative torque exerted by the other planet in a mean motion resonance. 

3.6. Analytic Considerations for e ~ 1 

We have seen how the timescales of the semimajor axis evolution and the eccentricity damping vary as we change 
the disk parameter. Especially, we have seen a strong dependence of t a and t e on the disk parameter a = p + (3 — q)/2 
if the planet eccentricity is high. In this section, we analytically investigate the behaviors of t a and t e in the limit of 
e -> 1. 

In this section, we approximate df/dt (see equation (|6ipl by 

± = «!m (74) 

dt r 2 ' [ ' 

since the contribution from the perturbing force is small compared to the force from the central star. We then obtain 



da_ 2(1 -e 2 ) 2 -" 
df ~ a 2 + a n 4 



e 2 sin 2 / + (1 + e cos ff' 2 { VI + e cos f - 1) 

(75) 



(1 + e cos /) 2 -« [e 2 sin 2 / + (1 + cos /)(2 + e cos / - 2^1 + e cos /)] 
and 

de H-e 2 ) 3 "" 1 

- -K 



3/2 



df a 3+a n 4 (l + e 2 + 2ecos/)(l + ecos/) 5 / 2 -« 

j sin 2 / [2e 2 (e + cos f)y/l + e cos / + e(l - e 2 )] 
\ [e 2 sin 2 /+ (1 + cos/)(2 + ecos/ - 2y/l + e cos /)] 3/2 
2(e + cos /) (1 + e cos f) 2 (y/2 + e cos / - 1) 1 



[e 2 sin 2 / + (1 + cos /) (2 + e cos / - 2Vl + ecos/)] 



3/2 



(76) 



From Figure [71 we see that the force exerted when the planet is near the aphelion and the perihelion is important 
in determining the orbital evolution of the planet. We therefore look at the places close to the perihelion or aphelion. 
Around these points, cos / and sin / are approximated by 

cos/~± (l + lsf) (77) 



2 

sin/~±<5/, (78) 

where we write f = Sf near the perihelion and / = ir + Sf near the aphelion, and the upper sign is for the perihelion 
and the lower sign is for the aphelion. We expand equations (|75]l and (f76|) up to the second order of Sf. After the 
expansion with respect to 8f, we take the limit of e — > 1. We define 

e = 1 - e (79) 

and take the lowest order of e. 
The calculations are tedious but straightforward. In the close vicinity of the perihelion, we obtain 

da e 2 - a 2(V2-1) 1 + 5f 2 (2 - V2)/8V2 

df ~ a 2 +«n 4 (3-2V2) 3 / 2 (l-<5/ 2 /4) 2 -«{l + (5/ 2 3(v / 2- l)/4(3-272)} 3 / 2 [ ' 

and _ _ _ 

de e 3 -" (V2-l)(l-<5/ 2 (5 + 6^/4^/2) m 

df~ a 3+Q n 4 V2(3-2V2)(l-<5/ 2 /4) 7 /2- Q (i + 5/2 3 (V2-l)/4(3-2y2)' 1 ' 
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We then integrate over the orbit. We denote the integral by, 

~da f 27V da 



df Jo df 

If we integrate over /, the terms in equations (fHUj) and (j8"Tj) that contain df result in a numerical factor of the order 
of unity. Therefore, close to the perihelion, da/ df and de/df averaged over one orbital period are 

^ oc a^ a e 2 - a (83) 



and 



§ oc a 3 - a e 3 ~ a . (84) 



We now turn the attention to the aphelion. In the close vicinity of the aphelion, we obtain 

da 1 1 - S f 2 If 3 / 2 

df a 2 + a n* (1 + Sf 2 /2e) 2 ~ a (l + Z5f 2 /2ef/ 2 



and 



de , „ 1 l + c5/ 2 /e 3 / 2 , x 

df a^n^s 1 / 2 (l + ^/ 2 /2e)5/2-«(l + ,5/2/ e 2)(i + 3 ( 5/2/2 £ )3/2- V ) 

We now integrate over / near the aphelion. The integration for da/df is rather complicated. We have 

f l-5f 2 /e 3 / 2 

Ia = J dSf (l + 5f 2 /2e) 2 -*{l + 35/ 2 /2e) 3 / 2 (87) 



and changing the integration variable to X = Sf/e 1 ^ 2 , 

la = e 1 ' 2 [ dX 



l-X 2 /e 1 ' 2 



(l + X 2 /2) 2 -«(l + 3X 2 /2) 3 / 2 ' 

If e is sufficiently small, the numerator is dominated by the second term, and therefore the resulting X a is negative 
and I a oc s . If e is not very small, the numerator is dominated by the first term. Therefore the resulting integral is 
positive and I a oc e 1 / 2 . Therefore, the contribution of the aphelion to da/df is 

da j a 4 ~ a negative for small e 
df \ a i ~ a e l l 2 positive for intermediate s 

The integration for de / df is more straightforward. We have 

j = f ,„ l+5f 2 /e 3 ' 2 

e J 1 {l + 5f 2 /2efl 2 - a {l+5f 2 /e 2 ){l + i5f 2 /2efl 2 { ' 

and changing the variable to X = Sf/e, 

/I _|_ £ l/2^2 
dX (1 + X 2 )(l + eX 2 /2f/ 2 - a {\ + ieX 2 /2f/ 2 ' (91) 

The terms proportional to e 1 / 2 and e are safely neglected, and therefore we have I e oc e. Therefore, for the contribution 
from the aphelion to de/df, we have 

% oc a 3 - a s^ 2 . (92) 
df 

It still remains to be determined whether contributions from the perihelion or from the aphelion dominate. The 
relative importance of the contributions from the aphelion and the perihelion depends on the disk model, and we 
numerically see which contribution is more important as we change a. In Figure [Til we plot the values of —da/df 
and —de/df for various disk models. We see that the contribution from the perihelion dominates when a > 2 
(corresponding to p ^ 1 if q = 1), and the expected behaviors of da/df oc e 2 ~ a and de/df oc e 3 ~ a are observed. If 

a <^ 2, the contribution from the aphelion comes into play, and we see de/df is roughly proportional to e 1 / 2 . We note 
that for the model with a — —I and a — 0, the direction of semimajor axis evolution is outward for the moderate values 
of eccentricity (e > 0.09 for a = model and e > 0.05 for a = —1), as we have already seen in the previous section. 
The inversion of the sign of da/df is also expected from the above analytic discussions at the aphelion. We have also 
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checked that the values of da/df are always positive if we have a large negative value of p, when the interaction is 
completely dominated by the contribution from the aphelion. 

We summarize the dependence of the semimajor axis evolution timescale and eccentricity damping timescale on the 
disk parameters in the case of highly eccentric planet. Since da/dt and da/df are related by 

da 1 da 



dt T orb df ' 



(93) 



t a oc a 5 / 2 da/df . In the same way, it is possible to show that t e oc ea 3 / 2 de/df oc a 3 l 2 dejdf , where we have used 
e ~ 1. In the case of a J> 2, the contribution from the perihelion dominates and therefore, 

t a oc a a - 3 / 2 e a - 2 (94) 
t e oc a a ~ 3/2 e a - 3 . (95) 

The evolution of the semimajor axis is always inward and the eccentricity always damps. If a <J 2, we take into account 
the contribution from the aphelion, 



, I a invvaiu lvl oiiicm c (96^ 

a 1 a a ~ 3 l 2 s~ x l 2 outward for intermediate e 

and 

t e oc a^^e- 1 ' 2 . (97) 

It is noted that the outward evolution of semimajor axis occurs only when a <^ and eccentricity always damps. The 
range of e where the planet experiences the outward semimajor axis evolution increases as the value of a is decreased. 
In the particular disk model with q = 1, a and p are related by a = p + 1. 

We have shown that there is a critical value for the power of the surface density, which determines the behavior of 
semimajor and eccentricity evolutions for a highly eccentric planet (e ~ 1). Let us consider the model with q = 1. If 
p > 2, both eccentricity and semimajor axis evolution timescale becomes small as we increase the planet's eccentricity. 
If 1 < p < 2, semimajor axis evolution timescale becomes shorter for higher eccentricity, while eccentricity evolution 
timescale becomes longer. For p < 1, both evolution timescales become longer as we increase the eccentricity of 
the planet. Further complication arises for the semimajor axis evolution for p <J 0, where the contribution from the 
aphelion comes into play to invert the direction of the evolution. 

However, we note that these critical values of p should not be overstated. In the model described above, the critical 
parameter that determines the behavior of the orbital evolution timescales is a = p+ (3 — q)/2. The appearance of this 
single parameter a is partly due to the prescription of the cutoff of the gravitational force. Note that the dynamical 
friction force we have used depends on the cutoff scale e and we have used the model in which e depends linearly on 
the disk scale height. The appearance of q in the parameter a depends on this simplified prescription of the cutoff. 
We al so note that the dependence of the timescales on the softening length is different from iPapaloizou and LarwoodI 
(200Cj), which is due to the difference in the formulation. 

Therefore, we consider that the critical values of p described above only have qualitative meanings. Yet, we expect 
that the model describes the qualitative behavior of the semimajor axis and eccentricity evolution of highly eccentric 
planets. The prediction is readily compared with two-dimensional numerical calculations, although a very large 
simulation box may be necessary. In order to get rid of the dependence on the softening parameter e, it is necessary 
to perform three-dimensional analyses. 

3.7. Validity of the Model 

In this section, we briefly discuss the applicability of our model. The most important assumption we have used 
in the model is the use of the supersonic dynamical friction formula (|41| . which is derived from linear perturbation 
analyses. The use of this formula assumes (1) the relative velocity is supersonic, (2) the flow in the vicinity of the 
planet is homogeneous, and (3) the steady state is reached instantaneously. We shall now discuss each assumption 
separately. We also discuss how non-linear effects on the dynamical friction force may change our results. 

The first assumption is justified if we consider a planet that has high enough eccentricity, as indicated in Figure OH 
Typically, when e 3> H/r, we can safely assume that the flow is supersonic. 

For the second assumption, let us consider which scale of the perturbation contributes to the force. In the discussion 
of Section [21 we have seen that the perturbation of all the scales larger than the cutoff length equally contributes to 
the dynamical friction force. Let us consider a planet embedded in a disk with Keplerian rotation, and the planet is 
located at the perihelion or aphelion. Then, there is a location in the disk where the relative velocity between the gas 
and the planet becomes zero. We call this radius an "instantaneous corotation radius", rc, and is given by, for the 
perihelion, 

GM* GM* 1 + e 

- "j ■ (9SV 

rc a 1 — e 

A similar equation can be derived for the aphelion. 



* 3 / 2 inward for small e 
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The assumption of the homogeneous flow breaks down if the length scale comparable to the distance between the 
planet's location and the instantaneous corotation radius is important in determining the force acting on the planet. If 
the planet's eccentricity is large, the distance between the instantaneous corotation radius rc and the planet's location 
is of the order of the scale of the disk itself (|rc — o\ ~ a ), which is much larger than the cutoff scale, which is comparable 
with the disk scale height. Therefore, in this case, we expect that the assumption of the homogeneous flow is justified. 
If the planet's eccentricity is small, the instantaneous corotation radius comes very close to the planet's orbital radius 
(rc ~ a), and the assumption of homogeneous medium becomes worse. Typically, we expect that this assumption 
breaks down if the planet's eccentricity is smaller than the disk aspect ratio, and in such a small eccentricity case, it 
is necessary to fully take into account the effects of shear around the planet's location. 

Regarding the third assumption, the time taken to reach the steady state may be estimated as follows. The length 
scale that is important in determining the force is of the order of the disk scale height, and the minimum speed of the 
propagation of the information of the perturbing potential may be the sound speed. This indicates that the time taken 
to develop the steady state is of the order of the Kepler timescale. However, since the relative velocity is supersonic, 
the time taken to develop the steady state is less than that, since the background flow will carry the information of 
the perturber. In any case, this assumption is only marginally satisfied. 

As discussed before, it is known from numerical simulations that the time when the maxima/mini ma of the torque 
occur s lags the time of the maxima/minima of the distance between the planet and the central star (e.g..fCrcsswcll ~t al.1 
2007). Our model does not give a precise prescription for this effects, and the instantaneous torque or power profiles 
are different from what we expect from the numerical simulations. However, if we take the average over one orbital 
period, our model is in good agreement with the previous calculations in the parameter range where both our model 
and the previous calculations are expected to give reasonable results. 

In short, we expect that our formulation is applicable to planets whose eccentricity is larger than the disk aspect 
ratio, in particular when one takes the average over one orbital period. 

Finally, we make a comment on the use of the dynamical friction formula derived by the linear perturbation analysis. 
iKim and K im (2009|) investigated the dynamical friction using the non-linear numerical simulations. They considered 
the case where the particle is embedded in a homogeneous, three-dimensional gas flow, and calculated the axisymmetric 
pattern of the gas flow induced by the particle's gravity. They found that there is a deviation from the linear theory 
depending on the Mach number Mq of the flow and the non-linear parameter 

A=^, (99) 
errs 

where rs is the (effective) radius of the body. The force exerted on the particle is deviated from the linear results by 
a factor of (77 /2) — °- 45 , where 77 = A/(Mq — 1). In the case of a planet embedded in a protoplanetary disk, A is of the 
order of 10 — 100, and therefore, non-linearity can be important. The correction factor, however, is the quantity of the 
order of unity, (since Mach number is of the order of 5-10) and therefore, our results may give a reasonable estimate 
of the timescale of the evolution of the orbital parameter. One possible effect of such non-linearity on our results is 
that the dependence of the timescales t a and t e on the orbital parameters of the planet can be different, since the force 
depends on the velocity of the perturbing potential in a different way. 

4. DISCUSSION: POSSIBLE APPLICATIONS TO OTHER STUDIES 

In this section, we discuss some possible applications of our results to other studies and the implications to the planet 
formation theory. We first show the fittin g formula for t a and t e derived from this study, which can be especially useful 
for population synthesis models such as llda and Linl f|2008T ) . We then discuss some implications from this study to 
recent planet formation scenarios. 

4.1. Fitting Formula for Timescale of Semimajor Axis and Eccentricity Evolution 
In this section, we derive the fitting formula for t a and t e . As seen before, within our framework, t a ,t e oc 
M~ 1 T,q 1 €oHo and the only disk parameter that controls the timescale is a = p + (3 — q)/2. Therefore, we fit 
the values of t a and t e as a function of the semimajor axis and the eccentricity for a given value of a with 1 < a < 4. 
For the data for fitting, we use the values obtained for 0.3 < a/ro < 10 and 0.1 < e < 0.99. 
Motivated by the results of the analytic considerations for e ~ 1, we use the fitting formula of the form: 

t „ . C M e<^ - W (^)" (±) (£) (100, 

Here, C a , Ca, A a and 77a are the function of a of the form 

C a (a) = C a0 ' (101) 

CaM = CaQ (£) ^ (102) 

a 



\ a (a) = X a0 [ - ) (103) 
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Va(a) = Vao ( — ) ■ (104) 



The same form of the fitting function is used for t e , 

U - C .(«,e- (1 - e-,-. (^)- (5^!)- (^) (£) , 



with 



• ' (^) " (106) 



Ce(o) = CeO ( - ) (107) 



a 



A e (o) = A e0 ( — ) (108) 



%(a) = Veo y-J ■ (109) 

We note that in the set of the fitting parameters, C a o and C e o have the dimension of time, while others are dimensionless. 

In Tables [1] and [21 we show the results of fitting for t a and t e , respectively. For C o and C e o, we give the values in 
the units of year with M* = M® and ro = 1AU. In the last column of the tables, we show the maximum values of the 
error of the fitting function, which is calculated as |t ca ic — ifltl Acalo where t ca i c is the value of t a or t e calculated in the 
model presented in this paper and t^t is the value of the fitting function. Our fitting formula successfully reproduce 
the values of t a and t e within 10 — 20% error. We also note that a Q_3 / 2 dependence that appears in equations (|M)) - (|97p 
can be clearly seen in the fitting parameters C ap and C ep . 

Here, we note that the caution must be payed in using the fitting formulae. In obtaining the fitting parameters, we 
have used the values with 0.3 < a/ro < 10 and 0.1 < e < 0.99. Especially, our formulation is not applicable to the case 
of a low-eccentricity planet (typically e <J H/r). It is necessary to combine the results of this paper and such formulae 
given by iTanaka and W ard ( 2004) which are applicable to the case of low eccentricity. 

It is also noted that the fitting formulae are given in the case of 1 < a < 4. Especially, our form of fitting functions 
(|100P and (|105[) fails to describe the change of the sign of t a that happens in the case of p <; —1 (see Figure ITTj) . 
Although we believe that this range of a covers most of the typical protoplanetary disk parameters, it is necessary to 
construct the fitting formula if one wants to apply for protoplanetary disk parameters out of the range of this fitting. 
Even in this case, however, we expect that the model described in this paper, which uses the dynamical friction force 
to calculate the evolution of orbital parameters, is still applicable. 

4.2. Implication to Planet Formation Theory 

We have developed a model for the disk-planet interaction that can be used especially for highly eccentric planets. 
We have found that in the case of a highly eccentric planet, the disk-planet interaction is essentially described by 
dynamical friction. We have derived how the evolution timescales of the planet's semimajor axis and eccentricity 
depend on disk parameters within our models. In this section, we discuss the possible implications to the planet 
formation t heory. 

Recently, iMachida et al.l ([20111 ) have suggested that a disk surrounding a newly born proto star is gravitationall y 
unstable, and therefore it may be possible to form a planet by gravitational instability (see also llnutsuka et aD feOlO). 
The planets born in such a way are expected to have high eccentricity. Also, it is possible to form a planet at a large 
distance from the central star. 

We briefly discuss the evolution timescale of such planets. From the fitting formula of the evolution timescales, it is 
possible to show that the orbital evolution timescales of the newly-born Jupiter-mass planets through the gravitational 
instability can be very rapid, since the disk is very massive compared to the central star. However, after the accretion 
phase onto the central star, the stellar mass becomes much larger than the disk mass. If there is a low-mass planet 
(e.g., sub- Jupiter mass) that has survived in the early phase and reasonable amount of eccentricity is m aintained, the 
orbit al evolution timescale can be of the comparable order of magnitude w ith the observed disk lifetime (Haisc h~et al.l 
2001). Such planets may be found by the direct imaging observations (e.g.. iMarois et all (|2008D ). although small mass 
planets may be very cold and the observations may be difficult. Nevertheless, a highly eccentric planet at a large orbital 
separation can be a good observational signature to test the theory of planet formation and disk-planet interaction. 

Another possible me c hanism to produce a planet i n a h ighly eccentric orbit is the planet-planet scattering (e.g., 
IChambers etHl (fl99fih :lM arzari and Weidenschilli ng (2002)). The scattering between planets can naturally result in 
the highly eccentric planets. From the calculations given in this paper, the evolution timescale of the orbital parameters 
can be very long even if the planet's orbital plane coincides with the gas disk plane. If the planet's orbital plane is 
not aligned with the plane of the disk, we expect that the evolution timescale becomes much larger, since the planet 
interacts with the disk only at a very small part of the orbit. The calculations given in this paper sets the lower limit 
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of the timescales of the orbital evolution for such planets. It is therefore possible that the scattered planets remain in 

a highly eccentric orbit, even if gas disk remains and the tidal interaction with the central star is not e ffective. 

An other interesting application of our model is the recently suggested "eccentricity trap" mechanism (Ogihara ct al. 
2010). They have suggested that if there is an inner cavity in a protoplanetary disk, it may be possible to maintain 
a planet at the inner edge of the disk. It is due to the positive torque acting on the planet at the aphelion, which 
stops the planet to migrate inward. Although our@formulation, which assumes a smooth disk profile, is not directly 
applicable to the disk with a sharp inner edge, we have already seen that the evolution of semimajor axis can be 
outward if there is a steep surface density gradient. It is an interesting extension to apply the method of dynamical 
friction to the disk with a cavity. 

5. SUMMARY 

In this paper, we have presented a model of the disk-planet interaction, which makes use of the dynamical friction 
force. This method can be especially useful for a highly eccentric planet. 

We have first derived the dynamical friction formula in a slab media using different methods from lOstrikerl (|1999f ). 
and calculated the dynamical friction force as a function of viscosity and Mach number. The integral that leads to the 
dynamical friction force is given by equation (|24[) , and some asymptotic expressions are given by equations (|28[) . (|30l) . 
(@DD, and 

We have then applied the dynamical friction formula to find the evolution of orbital parameters of an eccentric planet 
embedded in a disk. The results agree well with the previous calculations given bv iPapaloizou and Larwoodl (|2000H 
for moderate eccentricity (e ;> H/r), and we expect that it is possible to use our model to much higher eccentricity. 
In this sense, this paper provides a complement to the previous study. We have seen that the evolution of the orbital 
parameters depends on the disk structure, and calculated some critical values for p, the power of the surface density 
profile, at which the behaviors of the timescales of orbital evolutions are qualitatively altered. 

Our model presented in the paper is restricted to a planet with a high eccentricity, typically e ;> H/r, It is also 
necessary that the disk profile must be smooth in the vicinity of the planet, in order to justify the use of dynamical 
friction formula, which is obtained under the assumptions of a homogeneous slab. 
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APPENDIX 

TIME DEPENDENT ANALYSIS OF DYNAMICAL FRICTION FORCE 

In this appendix, we perform a time-dependent linear analysis for the gravitational interaction between a slab and 
an embedded particle. We do not perform the analysis in full details. The objective of this section is to show that it 
is sufficient to perform the steady-state analysis in order to obtain the dynamical friction force. 

We start with equation ([8|) , and perform Fourier transform in the spatial direction via equations (|16j) and f|lT[) . We 
then obtain an ordinary differential equation in time as 

d 2 a , 2 \ da 

We consider that at time t = 0, the particle is suddenly switched on, 



+ (2ik v va + vk 2 ) (c 2 fc 2 - k 2 y vl + ik 2 k y v v) a = -k 2 ip. (Al) 



(t < 0) 

- r 

-e- ke {t > 0) 



$(t,k) = { GM^_ ke ,^^ (A2) 



I k 

We define p, q, S (they should not be confused with p and q in the main text) as 

p = vk 2 + 2ik y VQ, (A3) 
q = c 2 fc 2 - vlk 2 + ik 2 k y v v, (A4) 

and 

S = -k 2 4> (A5) 

so equation (|Aip can be written in a concise way 

d 2 ct da . . 

_ +p _ +qa = S , (A6) 

where we have also removed the tilde sign for convenience. This equation should not be confused with that in the real 
space. 

If we define the function g by 

a{t)=e-v t ' 2 g{t), (A7) 
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equation (IA6|) becomes 



dt 2 



9 = S, 



(A8) 



where S = e pt / 2 S. 

The boundary condition is that there is no perturbation at t < 0. The solution that satisfies this condition is readily 
obtained as 

1 



9{t) 



g+(t) S(r)g-(r)dr-g-(t) 5(r) ff+ (r)dr 



where 



g± (*) = exp 

are the homogeneous solutions. From equation (|A9I) . we can write down the solution of a(t) as 

, . si 

a(t) = — < 1 ; = exp 



I { y/p 2 - 4q 

p/2-^/p 2 /4-q 



" f-VT-«l' 



\/p 2 - 4 <7 



exp 



P 
2 



P 



q\ t 



We first note that the arguments in the exponential in equation (|A11[) is given by 

1 



| ± -q = ^vk 2 + ik,,r t> i 



'^{kl + klf 



c 2 k 2 



(A9) 



(A10) 



(All) 



(A12) 



and therefore, the real part is always positive in the presence of viscosity. Therefore, in the limit of t — > oo, the 
contribution from the last two terms in equation (|A11[) vanishes if there is a viscosity, and the steady state solution 



S 



(A13) 



is reached. 

We now discuss the inviscid case, where the last two terms in equation (|A11[) oscillates. We show even in this case, 
the contribution to the dynamical friction force coming from these time-dependent terms vanishes as t . 

The expression of dynamical friction force (|20l) still holds in the time-dependent analysis. Therefore, we investigate 
the integral 

„-ke 



I, 



dkn. 



k 



■a. 



(A14) 



where the factor e ke jk comes from the Fourier transform of the gravitational potential. In the case of an inviscid 
gas, the surface density perturbation is given by 



*(*) = 



GMt 



c 2 k(l- M 2 k 2 /k 2 ) 
1 



x 1- 



M k v 



exp 



2"¥ lexp 



ick(l-M ^-)t 
k 



ick(l-M ^)t 
k 



+ exp 



exp 



-ickil + Mo-r)* 
k 



■ickil + M -Z)t 
k 



(A15) 



We consider the terms in the integral (|A14[) involving the time-dependence, which are second and the third terms in 
the curly bracket of equation (|A15I) . Changing the variables as done in equations (|22p and (1231) . we obtain the terms 
that are proportional to 



dk / d6 



sin" 9e- 2ke 
1 - M 2 sin 2 9 



exp [±ickt(l =F Afo sm ' 



(A16) 



where n = 1 for the second term and n = 2 for the third term. We further change the variable from k to ip = ckt, and 
perform the integration over ip. We now have 



Im(7 n ) = ±- 



d6 



sm 



ct J Q {(2e/ct) 2 + (1=f M sin6») 2 } (1 ± A/ sin 6») ' 



(A17) 
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where we have taken the imaginary part which contributes to the dynamical friction force. The integral involved in 
equation (|A17|) is finite for all the values of t > 0, and asymptotes to a finite value as t —> oo. Therefore, the contribution 
to the force arising from the time-dependent terms decays as t , which justifies the steady-state assumption given in 
the main text. 

THE EXACT EXPRESSION OF I„ 

The integral given by equation (|25|l can be performed analytically. We make use of Mathematica software to calculate 
the integral and we show only the results in this section. The integral is given by 

Ie(p,q) = ^, (Bl) 



where the numerator is 
A 
2 



4 = vV (-VvVl - 4P -2p + q + 2 



x log 



2p 2 - - Ap -2p + q 



'2p 3 + p 2 (y/qy/q - Ap - q - 2) + q 3 l 2 yfq - Ap + p (Aq - 2^/qy/q - Ap) - q 2 ( 
+V2p 2 (y/qy/q-4p + 2p-q-2) 



x log 



-2p 2 + y/qy/q^Ap + 2p - q 



2p 3 + P 2 (y/qVl ~ 4 P - 9 - 2) + q 3 / 2 x/q-Ap + p (Aq - 2^/qy/q - Ap) - q 2 



+TT[VqVq-Ap-2p + q) (B2) 



xJ2p 3 +p 2 ( y /q~y/q-Ap-q-2) + q 3 / 2 yj q - Ap + p [Aq - 2y/q~ yj q - Ap] -q 



,2 



x / -2p 4 + 4p 3 - 2p 2 (yJqy/q~A3 + q + l) - q 3 / 2 (V^"4^ + yjq) + 2p {^qVq~~A^ + 2q) 
X y 2p 3 - P *( y /qy/q-=^ + q + 2) - q 3 / 2 {Vq=4p + y/q) + 2p + 2?) 

and the denominator is 

# = VV? - 4 P ( 2 P 2 + - 2p + q) 

X ^2p 3 + p 2 (y/q-yjq - Ap -q-2) 

\ 1/2 



+q 3 / 2 y/q~^Ap~ + p (iq - 2^^q~Ap^j - q 2 . (B3) 
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Fig. 2. — The dynamical friction force for Re = 100. The results of numerical integration of equat ion 12 4 1) are shown by plus symbols, 
and the analytic formulae in the limit of supersonic and subsonic cases, which are given by equations J30D and (14 1 I t respectively, are shown 
by dashed lines. 
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Fig. 3. — Definition of unit vectors ey and ejv and their relations with the unit vectors in the cylindrical coordinate e r and e^. The 
angle between e r and is j3. 
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Fig. 4. — The evolution of relative velocity vector along the planetary orbit with e = 0.5. Solid line shows the orbit of the planet, and 
the arrows show the direction of relative velocity vector between the gas and the planet. We have assumed that the central star is at the 
origin, and the planet and the disk are both rotating in the counterclockwise direction. 
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Fig. 6. — The power V (left panel) and torque T (right panel) exerted on the planet with eccentricity e = 0.3 (solid line) and e = 0.5 
(dashed line). The fiducial model is used. For the power, we plot (V /nE p ) / (M p / M*) and for the torque, we plot (T/nL p )/(M p /M*), 
where E p and L p are the planet's energy and angular momentum, respectively, and n is the mean motion of the planet. The horizontal 
axis shows the true anomaly /, and / = corresponds to the perihelion. 
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Fig. 7. — The evolution of semimajor axis da/dt (left panel) and eccentricity de/dt (right panel) of the planet with eccentricity e = 0.3 
(solid line) and e = 0.5 (dashed line). The fiducial model is used. For the evolution of semimajor axis, we plot (da/dt)(l/na) / (Mp /M*) 
and for the evolution of eccentricity, we plot (de/di)(l/n)/ (M p /M»). The horizontal axis shows the true anomaly /, and / = corresponds 
to the perihelion. 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Fig. 8. — The migration timescale t m (left) and the eccentricity damping timescale t e (right) for fiducial model. For both of the panels, 
we compare our results (crosse s) and the formula by IFapaloizou and Larwoodl J2000I ) (solid jine). For t m , we al s o sho w the formula of 
Papaloizou and Larwood (2000) multiplied by 3/2 is shown by dashed line for e < 0.5, which Cresswell and Nelson (2006) claims to fit the 
results of numerical simulations better. 
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Fig. 9. — The timescale of the evolution of semimajor axis t a (left) and the timescale of the evolution of eccentricity t e (right) for the 
models with a = 1,2,3,4 (p = 0,1,2,3, respectively for q = 1). The lines for a = 1 and a = 2 are almost identical for t e . We take 
a / r = 1 an< i the unit of the time is year with ro = 1AU and M t = Mq. The direction of the semimajor axis evolution is inward, and the 
eccentricity always damps. 
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Fig. 10. — The timescale of the evolution of semimajor axis t a (left) and the timescale of the evolution of eccentricity t e (right) for the 
models with a = —1,0, 1 (p = —2, —1,0, respectively for q = 1). We take a/ro = 1 and the unit of the time is year with ro = 1AU and 
M* = Mq. The left panel shows l/t a for clarity, since the sign of the direction of the semimajor axis evolution changes. The negative sign 
indicates that the direction of the semimajor axis evolution is outward. Eccentricity always decrease for the parameters presented here. 
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Fig. 11. — The orbit-averaged values of —da/df (left) and —de/df (right) for various eccentricity and disk parameters with a/ro = 1- 
The power p (the parameter a = p + 1 in the case of q = 1) of the surface density is varied and denoted in the figures. The horizontal axis 
shows the values of 1 — e. For da/df, the absolute values are plotted if —da/df is negative. Only the cases when a = —1 and a = with 

large values of 1 — e (1 — e > 0.09 for a = and 1 — e > 0.05 for a = —1), wc obtain negative values of —da/df . For the values of —de/df, 
the values are always positive. 
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TABLE 1 
Fitting results for t a 



a 


C a0 [yr] 


Cap 


CaO 


Cap 


AaO 


^ap 


VaO 




max. error 


1.0 


3.789E+06 


-5.241E-01 


9.759E-01 


-1.471E-05 


1.320E-03 


-1.211E-01 


-1.986E-01 


3.019E-09 


1.886E-01 


1.5 


3.762E+06 


4.353E-07 


8.977E-01 


3.777E-08 


1.658E-01 


2.010E-06 


-2.173E-01 


1.068E-09 


3.895E-02 


2.0 


3.988E+06 


5.000E-01 


1.008E+00 


-1.188E-09 


1.129E+00 


7.216E-08 


5.318E-02 


3.436E-09 


4.497E-03 


2.5 


2.861E+06 


1.000E+00 


9.952E-01 


0.000E+00 


2.112E+00 


0.000E+00 


4.984E-01 


0.000E+00 


2.038E-03 


3.0 


2.265E+06 


1.500E+00 


9.944E-01 


0.000E+00 


2.063E+00 


0.000E+00 


9.957E-01 


0.000E+00 


2.787E-03 


3.5 


1.948E+06 


2.000E+00 


1.010E+00 


1.604E-10 


1.920E+00 


-3.424E-10 


1.503E+00 


-4.489E-11 


3.085E-03 


4.0 


1.793E+06 


2.500E+00 


1.041E+00 


6.360E-10 


1.774E+00 


-1.094E-09 


2.012E+00 


-1.386E-10 


1.058E-02 



TABLE 2 
Fitting results for t e 



a 


C e o [yr] 


Cep 


CeO 


Cep 


A e 


Aep 


VeO 




max. error 


1.0 


3.426E+04 


-5.090E-01 


2.651E+00 


-1.908E-06 


1.158E-03 


-1.331E-02 


-6.765E-01 


8.473E-10 


8.247E-02 


1.5 


3.768E+04 


-9.810E-02 


2.615E+00 


-3.975E-05 


2.431E-03 


-1.324E-01 


-7.394E-01 


-2.246E-08 


1.126E-01 


2.0 


3.307E+06 


5.000E-01 


2.914E+00 


0.000E+00 


1.258E+00 


-9.885E-11 


-7.597E-01 


0.000E+00 


5.556E-02 


2.5 


3.956E+06 


1.000E+00 


2.984E+00 


0.000E+00 


1.733E+00 


-4.559E-11 


-5.236E-01 


0.000E+00 


9.257E-03 


3.0 


2.372E+06 


1.500E+00 


2.991E+00 


0.000E+00 


-7.393E-09 


-5.798E-09 


-2.964E-02 


0.000E+00 


1.802E-02 


3.5 


4.528E+06 


2.000E+00 


3.038E+00 


0.000E+00 


1.376E+00 


-5.197E-10 


5.102E-01 


0.000E+00 


9.239E-03 


4.0 


7.059E+06 


2.500E+00 


3.183E+00 


0.000E+00 


8.549E-01 


-5.166E-10 


1.027E+00 


0.000E+00 


2.573E-02 



